function [beta_new,r] = newton_raphson(epsNR,maxNR,n,m,G,Gprime,Gprimeprime,...
        weight_function,sup_norm,beta_old,y,v)
    
% Function which performs the Newton-Raphson iteration to find beta
% Accepts either weighting function (good or bad)

% Initialize
%
p(1:n) = NaN;
deriv_sum = zeros(m,m);
r=1;
beta_diff = 1e10;

% Newton-Raphson iteration
%
% Perform while number of iterations is less than a specified value
% (maxNR), and the difference in the calculated betas from the previous two
% steps is greater than a specified value (epsNR)
while r<=maxNR && beta_diff>epsNR
    vector_sum = zeros(m,1);
    
    % Compute the vector term
    for i = 1:n
        p(i) = G(beta_old * v(i,:)');
        vector_sum = vector_sum + Gprime(p(i)) * v(i,:)' ...
            * weight_function(p(i),y(i));
    end
    
    % Compute the matrix (derivative) term
    for kk = 1: m
        for jj = 1: m
            deriv_sum(kk,jj) = 0;
            for i = 1: n
                wp = weight_function(p(i),y(i));
                Gp = Gprime(p(i));
                deriv_sum(kk,jj) = deriv_sum(kk,jj) + ...
                    (wp * Gprimeprime(p(i)) - wp^2 * Gp) ...
                    * Gp * v(i,kk)*v(i,jj);
            end
        end
    end
    
    % Invert matrix
    inverse_deriv = inv(deriv_sum);
    
    % Newton-Raphson step
    beta_new = beta_old - transpose(inverse_deriv * vector_sum);
    
    % Stop if already diverged
    if(sum(isnan(beta_new))>0)
        r = NaN;
    end
    
    % Prepare for next step
    beta_diff = sup_norm(beta_new,beta_old);
    beta_old = beta_new;
    r = r+1;
end
